Analysis of The National Health and Nutrition Examination Survey (NHANES)

1 INTRODUCTION

nhanes <- read.csv("https://jmusa3.github.io/jmusa/nhanes.csv")

#CREATE MISSING VALUES in some variables

# Set seed for reproducibility
set.seed(123)  

# Determine the number of missing values (15% of total rows)
num_missing <- round(0.14 * nrow(nhanes))

# Randomly select row indices to introduce missing values 
missing_indices <- sample(1:nrow(nhanes), num_missing, replace = FALSE)

# Set selected rows in body.weight to NA
nhanes$body.weight[missing_indices] <- NA

# Verify missing values
#sum(is.na(nhanes$body.weight))  


# Set seed for reproducibility
set.seed(456)  

# Determine the number of missing values (15% of total rows)
num_missing <- round(0.15 * nrow(nhanes))

# Randomly select row indices to introduce missing values
missing_indices <- sample(1:nrow(nhanes), num_missing, replace = FALSE)

# Set selected rows in body.weight to NA
nhanes$cholestrol[missing_indices] <- NA


# Set seed for reproducibility
set.seed(789)  

# Determine the number of missing values (15% of total rows)
num_missing <- round(0.15 * nrow(nhanes))

# Randomly select row indices to introduce missing values
missing_indices <- sample(1:nrow(nhanes), num_missing, replace = FALSE)

# Set selected rows in race to NA
nhanes$race[missing_indices] <- NA

# Verify missing values

#sum(is.na(nhanes$race))  
#sum(is.na(nhanes$cholestrol)) 
#sum(is.na(nhanes$body.weight)) 

# create a copy of nhanes for MICE
nhanes1 <- nhanes
nhanes1$race <- as.factor(nhanes1$race)
nhanes1$cholestrol <- as.numeric(nhanes1$cholestrol)

The National Health and Nutrition Examination Survey (NHANES) is a program conducted by the National Center for Health Statistics (NCHS), a division of the Centers for Disease Control and Prevention (CDC). NHANES has been conducted in various forms since the early 1960s, with the current continuous survey format beginning in 1999. The survey is designed to assess the health and nutritional status of adults and children in the United States through a combination of interviews, physical examinations, and laboratory tests. Data is collected from a nationally representative sample of the U.S. population using a complex, multistage probability sampling design. Participants undergo household interviews followed by health examinations at mobile examination centers (MECs), where trained medical professionals collect biological samples and conduct clinical assessments. The NHANES data provides critical insights into public health trends, helping policymakers, researchers, and healthcare professionals make informed decisions on national health priorities. The purpose of this study is to assess the health and nutritional status of the study participants with metrics such as high blood pressure.

There are 15 features and 7927 observations in this data set. These features are listed below.

1. obs: Respondent Identification Number
2. psu: Pseudo-PSU 1,2 psu
3. stratum: Pseudo-stratum (01 - 49)
4. stat.weight: Statistical weight (225.93 - 139744.9)
5. age: Age (years)
6. sex: 0 = Female, 1 = Male sex
7. race: Race (1 = White; 2 = Black; 3 = Other)
8. body.weight: Body Weight (pounds)
9. height: Standing Height (inches)
10. avg.sbp: Average Systolic BP (mm/Hg)
11. avg.dbp: Average Diastolic BP (mm/Hg)
12. past.smk: Has respondent smoked > 100 cigarettes in life (1 = Yes, 2 = No)
13. current.smk: Does respondent smoke cigarettes now? (1 = Yes, 2 = No)
14. cholestrol: Serum Cholesterol (mg/100ml)
15. hbp: High Blood Pressure (0 if PEPMNK1R <= 140; 1 if PEPMNK1R > 140)

The obs, psu, stratum, and stat.weight variables are not included in this analysis because they do not directly provide health-related information.

The data set has 317 missing values for body weight and 396 missing values for cholesterol.

2 EDA

2.1 Distribution of Individual Features

The distributions of the individual features in the NHANES dataset provide a snapshot of the characteristics of the population in terms of various health and demographic variables. Understanding these distributions is essential for identifying patterns, detecting potential biases, and assessing the range and variability of each feature. This analysis helps in understanding how variables such as age, body weight, blood pressure, and lifestyle factors are spread across the dataset, which is crucial for subsequent statistical modeling and interpretation.

# Convert variables to factors with proper labels
nhanes$sex <- factor(nhanes$sex, levels = c(0, 1), labels = c("Female", "Male"))
nhanes$race <- factor(nhanes$race, levels = c(1, 2, 3), labels = c("White", "Black", "Other"))
nhanes$past.smk <- factor(nhanes$past.smk, levels = c(1, 2), labels = c("Yes", "No"))
nhanes$current.smk <- factor(nhanes$current.smk, levels = c(1, 2), labels = c("Yes", "No"))

# Create frequency tables
sex_count <- table(nhanes$sex)

race_count <- table(nhanes$race)

past.smk_count <- table(nhanes$past.smk)

current.smk_count <- table(nhanes$current.smk)


par(mfrow = c(2, 2))

# Plot the bar charts
barplot(height = sex_count, col = "steelblue", main = "Sex", xlab = "Sex", ylab = "Count")

barplot(height = race_count, col = "steelblue", main = "Race", xlab = "Race", ylab = "Count")

barplot(height = past.smk_count, col = "steelblue", main = "Past Smoker", xlab = "Past Smoker", ylab = "Count")

barplot(height = current.smk_count, col = "steelblue", main = "Current Smoker", xlab = "Current Smoker", ylab = "Count")

There are slightly more male observations than female observation in this data set. The majority of participants are white, followed by black. Most participants were past smokers. In fact, there is an imbalance for past smoker. About half of current participants are current smokers.

par(mfrow = c(2, 2))  

# Histogram for Cholesterol
hist(nhanes$cholestrol, 
     main = "Histogram of Cholesterol", 
     xlab = "Cholesterol", 
     col = "steelblue", 
     border = "black", 
     breaks = 20)  # Control number of bins
    
     
# Histogram for Systolic BP
hist(nhanes$avg.sbp, 
     main = "Histogram of Systolic BP", 
     xlab = "Systolic BP", 
     col = "steelblue", 
     border = "black", 
     breaks = 20)
    

# Histogram for Body Weight
hist(nhanes$body.weight, 
     main = "Histogram of Body Weight", 
     xlab = "Body Weight", 
     col = "steelblue", 
     border = "black", 
     breaks = 20)
  
     
# Histogram for Height
hist(nhanes$body.weight, 
     main = "Histogram of Height", 
     xlab = "Height", 
     col = "steelblue", 
     border = "black", 
     breaks = 20)   

Most participants do not have high blood pressure. The age of participants varies from about 20 to about 90. The height distribution is roughly symmetric. The body weight distribution is skewed right due to some high outliers. The high outliers may indicate obesity.

par(mfrow = c(2,2))

hist(nhanes$avg.sbp, col = "steelblue", 
     main = "Systolic BP",
     xlab = "Systolic BP", ylab = "Frequency",
     breaks = 30)

hist(nhanes$avg.dbp, col = "steelblue", 
     main = "Diastolic BP",
     xlab = "Diastolic BP", ylab = "Frequency",
     breaks = 30)

hist(nhanes$cholestrol, col = "steelblue", 
     main = "Cholesterol",
     xlab = "Cholesterol", ylab = "Frequency",
     breaks = 30)

Systolic BP and cholesterol have skewed right distributions. Diastolic BP is roughly symmetric. Of the 7927 participants, 1964 are obese.

2.2 Relationship between Features

Examining relationships between variables using scatterplots and a correlation matrix provides a comprehensive approach to understanding associations. Scatterplots visually reveal how two continuous variables interact, allowing us to identify trends, correlations, and potential outliers. Meanwhile, a correlation matrix quantifies the strength and direction of relationships between multiple variables at once, helping to pinpoint which pairs are strongly correlated. Together, scatterplots and a correlation matrix offer valuable insights into the underlying patterns in the data, guiding further statistical analysis and model selection.

# Set seed for reproducibility
set.seed(123)

# Sample 500 rows (adjust the number as needed)
sample_indices <- sample(nrow(nhanes), size = 500, replace = TRUE)
sampled_data <- nhanes[sample_indices, c("age", "body.weight", "height", "avg.sbp", "avg.dbp", "cholestrol")]

# Create a pairwise scatterplot matrix with the sampled data
pairs(sampled_data, 
      main = "Pairwise Scatterplots (Bootstrapped Sample)",
      pch = 19,           # Shape of the points (circle)
      col = rgb(0, 0, 1, 0.5)  # Color with transparency (blue)
)

nhanes_subset <- nhanes[, c("age", "body.weight", "height", "avg.sbp", "avg.dbp", "cholestrol")]

# Compute the correlation matrix
cor_matrix <- cor(nhanes_subset, use = "pairwise.complete.obs")

# Print the correlation matrix
print(cor_matrix)

                    age body.weight      height    avg.sbp   avg.dbp cholestrol
age          1.00000000 -0.01377966 -0.08538971 0.55641362 0.1143690  0.2738272
body.weight -0.01377966  1.00000000  0.45314709 0.13693751 0.2900425  0.1032280
height      -0.08538971  0.45314709  1.00000000 0.01379185 0.2060106 -0.0631896
avg.sbp      0.55641362  0.13693751  0.01379185 1.00000000 0.5213667  0.2306877
avg.dbp      0.11436905  0.29004249  0.20601064 0.52136671 1.0000000  0.1706567
cholestrol   0.27382718  0.10322800 -0.06318960 0.23068766 0.1706567  1.0000000

corrplot(cor_matrix, method = "color", type = "upper", tl.col = "black", tl.cex = 0.8)

The most highly correlate features are age and average systolic blood pressure, body weight and height, average systolic blood pressure and average diastolic blood pressure, age and cholesterol, body weight and cholesterol, and body weigh and average diastolic blood pressure. The pairwise relationships are displayed with individual scatter plots.

# Set seed for reproducibility
set.seed(123)

# Bootstrap a smaller sample (adjust size as needed)
sample_indices <- sample(nrow(nhanes), size = 500, replace = TRUE)
sampled_data <- nhanes[sample_indices, ]

# Set up a 2x3 plotting grid
par(mfrow = c(2, 3))

# Scatterplot for Age and Avg. Systolic Blood Pressure
plot(sampled_data$age, sampled_data$avg.sbp, main = "Age vs Avg. Systolic BP", 
     xlab = "Age", ylab = "Avg. Systolic BP", col = "blue", pch = 19)

# Scatterplot for Body Weight and Height
plot(sampled_data$body.weight, sampled_data$height, main = "Body Weight vs Height", 
     xlab = "Body Weight", ylab = "Height", col = "red", pch = 19)

# Scatterplot for Avg. Systolic Blood Pressure and Avg. Diastolic Blood Pressure
plot(sampled_data$avg.sbp, sampled_data$avg.dbp, main = "Avg. Systolic vs Avg. Diastolic BP", 
     xlab = "Avg. Systolic BP", ylab = "Avg. Diastolic BP", col = "forestgreen", pch = 19)

# Scatterplot for Age and Cholesterol
plot(sampled_data$age, sampled_data$cholestrol, main = "Age vs Cholesterol", 
     xlab = "Age", ylab = "Cholesterol", col = "purple", pch = 19)

# Scatterplot for Avg. Systolic BP and Cholesterol
plot(sampled_data$avg.sbp, sampled_data$cholestrol, main = "Avg Systolic BP vs Cholesterol", 
     xlab = "Cholesterol", ylab = "Avg.SBP", col = "gold", pch = 19)

# Scatterplot for Body Weight and Avg. Diastolic Blood Pressure
plot(sampled_data$body.weight, sampled_data$avg.dbp, main = "Body Weight vs Avg. Diastolic BP", 
     xlab = "Body Weight", ylab = "Avg. Diastolic BP", col = "maroon", pch = 19)

Six of the most correlated pairwise relationship are visualized in scatter plots. All pairwise relationships are linear with no apparent clustering.

# Set up a 2x2 plotting grid
par(mfrow = c(2, 2))

# Create a contingency table
sex_hbp_table <- table(nhanes$sex, nhanes$hbp)

# Generate mosaic plot
mosaicplot(sex_hbp_table, 
           col = c("gold", "royalblue"), 
           main = "Sex and High Blood Pressure",
           xlab = "Sex Status", 
           ylab = "HBP Status",
           border = "black")


# Create a contingency table
race_hbp_table <- table(nhanes$race, nhanes$hbp)

# Generate mosaic plot
mosaicplot(race_hbp_table, 
           col = c("gold", "royalblue"), 
           main = "Race and High Blood Pressure",
           xlab = "Race", 
           ylab = "HBP Status",
           border = "black")


# Create a contingency table
current.smk_hbp_table <- table(nhanes$current.smk, nhanes$hbp)

# Generate mosaic plot
mosaicplot(current.smk_hbp_table, 
           col = c("gold", "royalblue"), 
           main = "Current Smoking Status and HBP",
           xlab = "Current Smoking Status", 
           ylab = "HBP Status",
           border = "black")


# Create a contingency table
past.smk_hbp_table <- table(nhanes$past.smk, nhanes$hbp)

# Generate mosaic plot
mosaicplot(past.smk_hbp_table, 
           col = c("gold", "royalblue"), 
           main = "Past Smoking Status and HBP",
           xlab = "Past Smoking Status", 
           ylab = "HBP Status",
           border = "black")

Relationship between high blood pressure and race, sex, and current smoking status can be visualized in the mosaic plots. There appears to be an association between high current smoking and race with high blood pressure. A less pronounced relationship can been seen between sex and race with high blood pressure.

# Set up a 2x2 plotting grid
par(mfrow = c(1, 2))


# Create a contingency table
current.smk_race_table <- table(nhanes$race, nhanes$current.smk)

# Generate mosaic plot
mosaicplot(current.smk_race_table, 
           col = c("gold", "royalblue"), 
           main = "Race and Current Smoke Status",
           xlab = "Race", 
           ylab = "Current Smoke",
           border = "black")


# Create a contingency table
current.smk_sex_table <- table(nhanes$sex, nhanes$current.smk)

# Generate mosaic plot
mosaicplot(current.smk_sex_table, 
           col = c("gold", "royalblue"), 
           main = "Sex and Current Smoke Status",
           xlab = "Sex", 
           ylab = "Current Smoke",
           border = "black")

# Create a contingency table
past.smk_race_table <- table(nhanes$race, nhanes$past.smk)

An association exists between race and current smoking. Black individuals tend to smoke more than white or other individuals. An association exists between sex and current smoke status. Female tend to smoke currently slight more than males in this population.

3 MISSING VALUES

3.1 Imputation for Categorical Features

The distribution of missing values for is examined to determine any patterns. This step is done to ensure that value are missing at random.

# Set up a 2x2 plotting grid
par(mfrow = c(2, 2))

#sum(is.na(nhanes$race))

#mean(is.na(nhanes$race))

#table(nhanes$race, useNA = "ifany")


#ggplot(nhanes, aes(x = factor(race))) +
  #geom_bar(fill = "blue") +
  #labs(title = "Distribution of Race (Before Imputation)", x = "Race", y = "Count") +
  #theme_minimal()


aggr(nhanes, col = c("navyblue", "red"), numbers = TRUE, sortVars = TRUE, 
     labels = names(nhanes), cex.axis = 0.7, gap = 2, 
     ylab = c("Missing Data Overview", "Pattern of Missing Data"))

 Variables sorted by number of missings: 
    Variable     Count
        race 0.1499937
  cholestrol 0.1499937
 body.weight 0.1400278
           X 0.0000000
         obs 0.0000000
         psu 0.0000000
     stratum 0.0000000
 stat.weight 0.0000000
         age 0.0000000
         sex 0.0000000
      height 0.0000000
     avg.sbp 0.0000000
     avg.dbp 0.0000000
    past.smk 0.0000000
 current.smk 0.0000000
         hbp 0.0000000

A pattern of missing values exists between body weight, BMI, and obese. This pattern can be explained by the fact that body weight is used to calculate BMI, which is then used to classify observations for obese. However, the pattern of missing value appears to be random for the other features.

No pattern of missing values exists between race and other features, In turn, kNN imputation is used for missing values of race.

sum(is.na(nhanes$race))

[1] 1189

table(nhanes$race, useNA = "ifany")  # Check distribution including NAs

White Black Other  <NA> 
 4693  1858   187  1189 

# Perform KNN Imputation (single imputation)
nhanes_imputed <- kNN(nhanes, variable = "race", k = 3)

# Keep only the original columns (without extra columns added by kNN)
nhanes_imputed <- nhanes_imputed[, names(nhanes)]



# COMPARE before and after imputation distribution

par(mfrow = c(1,2))  # Split plotting area into two

# Convert race to a factor for better visualization
nhanes$race <- as.factor(nhanes$race)
nhanes_imputed$race <- as.factor(nhanes_imputed$race)

# Create a new column to indicate whether the data is original or imputed
nhanes$source <- "Before Imputation"
nhanes_imputed$source <- "After Imputation"

# Combine both datasets
nhanes_compare <- bind_rows(nhanes, nhanes_imputed)

# Plot the distribution before and after imputation
ggplot(nhanes_compare, aes(x = race, fill = source)) +
  geom_bar(position = "dodge") +
  labs(title = "Distribution of Race Before and After Imputation",
       x = "Race",
       y = "Count",
       fill = "Dataset") +
  theme_minimal()

imputed_nhanes <- nhanes_imputed

The overall distribution of race after imputation changed slightly after imputation.

3.2 Regression-Based Imputation for Numerical Feature

The pattern of missing values for body weight and cholesterol is examined for any potential patterns.

#vis_miss(nhanes[, c("body.weight", "cholestrol")])


# Set smaller font size for the plot
par(cex = 0.4)  # Adjust this value for desired font size

# Generate the missing data pattern plot

md.pattern(imputed_nhanes)

     X obs psu stratum stat.weight age sex race height avg.sbp avg.dbp past.smk
5792 1   1   1       1           1   1   1    1      1       1       1        1
1025 1   1   1       1           1   1   1    1      1       1       1        1
946  1   1   1       1           1   1   1    1      1       1       1        1
164  1   1   1       1           1   1   1    1      1       1       1        1
     0   0   0       0           0   0   0    0      0       0       0        0
     current.smk hbp source body.weight cholestrol     
5792           1   1      1           1          1    0
1025           1   1      1           1          0    1
946            1   1      1           0          1    1
164            1   1      1           0          0    2
               0   0      0        1110       1189 2299

# Reset graphics settings (optional)
par(cex = 1)

A pattern of missing values exists between body weight, BMI, and obese because

Body weight is correlated with height, average diastolic blood pressure, and sex. In turn, the missing body weight values will be imputed using a regression model using said explanatory features.

# Subset nhanes to contain only complete observations for selected variables
complete_nhanes <- imputed_nhanes[complete.cases(imputed_nhanes[, c("body.weight", "height", "avg.dbp", "sex")]), ]

# Fit the linear regression model
model_body.weight <- lm(body.weight ~ height + avg.dbp + sex, data = complete_nhanes)

# View the model summary
summary(model_body.weight)

Call:
lm(formula = body.weight ~ height + avg.dbp + sex, data = complete_nhanes)

Residuals:
    Min      1Q  Median      3Q     Max 
-86.131 -22.276  -4.259  17.073 269.863 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -207.90214    9.27828 -22.407  < 2e-16 ***
height         4.86582    0.14277  34.080  < 2e-16 ***
avg.dbp        0.76728    0.03913  19.606  < 2e-16 ***
sexMale       -7.00000    1.07917  -6.486  9.4e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 33.52 on 6813 degrees of freedom
Multiple R-squared:  0.2495,    Adjusted R-squared:  0.2491 
F-statistic: 754.9 on 3 and 6813 DF,  p-value: < 2.2e-16

# Plot the residuals
plot(model_body.weight$residuals)

The residual plot shows no patterns. In turn, the model fit the data well and the assumption of constant variance of residuals is satisfied.

# Create a copy of the dataset to work with
#imputed_nhanes <- complete_nhanes

colSums(is.na(imputed_nhanes))

          X         obs         psu     stratum stat.weight         age 
          0           0           0           0           0           0 
        sex        race body.weight      height     avg.sbp     avg.dbp 
          0           0        1110           0           0           0 
   past.smk current.smk  cholestrol         hbp      source 
          0           0        1189           0           0 

# Identify the rows with missing values for body.weight
missing_body_weight <- is.na(imputed_nhanes$body.weight)

# Use the model to predict the missing values
imputed_nhanes$body.weight[missing_body_weight] <- predict(model_body.weight, 
                                                           newdata = imputed_nhanes[missing_body_weight, c("height", "avg.dbp", "sex", "age")])

# Check the dataset
#head(imputed_nhanes)

# Create a plot to overlay the density curves
ggplot() +
  geom_density(data = complete_nhanes, aes(x = body.weight), fill = "blue", alpha = 0.2, na.rm = TRUE) + 
  geom_density(data = imputed_nhanes, aes(x = body.weight), fill = "red", alpha = 0.2, na.rm = TRUE) +
  labs(title = "Density Plot of Body Weight: Before and After Imputation",
       x = "Body Weight", y = "Density") +
  theme_minimal() +
  scale_fill_manual(name = "Body Weight", values = c("blue" = "blue", "red" = "red"), 
                    labels = c("Before Imputation", "After Imputation"))

The distribution before and after imputation is similar. In turn, the imputation did not change the distribution of body weight.

The features that are most correlated to cholesterol are age, average systolic blood pressure, and body weight. In turn, the missing values of cholesterol will be imputed with a model using these features.

# Subset nhanes to include only rows with complete cases for cholestrol, age, avg.sbp, and BMI
complete_nhanes <- imputed_nhanes[!is.na(imputed_nhanes$cholestrol) & !is.na(imputed_nhanes$age) & 
                           !is.na(imputed_nhanes$avg.sbp) & !is.na(imputed_nhanes$body.weight), ]

complete_nhanes <- imputed_nhanes[complete.cases(imputed_nhanes[, c("cholestrol", "age", "avg.sbp", "body.weight")]), ]

# Fit the linear regression model to predict cholestrol
model_cholestrol <- lm(cholestrol ~ age + avg.sbp + body.weight, data = complete_nhanes)
summary(model_cholestrol)

Call:
lm(formula = cholestrol ~ age + avg.sbp + body.weight, data = complete_nhanes)

Residuals:
    Min      1Q  Median      3Q     Max 
-162.16  -28.18   -3.08   24.08  498.54 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 135.10667    3.96327  34.090  < 2e-16 ***
age           0.54613    0.03465  15.762  < 2e-16 ***
avg.sbp       0.21908    0.03267   6.707 2.15e-11 ***
body.weight   0.10341    0.01437   7.198 6.76e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 42.59 on 6734 degrees of freedom
Multiple R-squared:  0.09077,   Adjusted R-squared:  0.09037 
F-statistic: 224.1 on 3 and 6734 DF,  p-value: < 2.2e-16

# View the model summary
summary(model_cholestrol)

Call:
lm(formula = cholestrol ~ age + avg.sbp + body.weight, data = complete_nhanes)

Residuals:
    Min      1Q  Median      3Q     Max 
-162.16  -28.18   -3.08   24.08  498.54 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 135.10667    3.96327  34.090  < 2e-16 ***
age           0.54613    0.03465  15.762  < 2e-16 ***
avg.sbp       0.21908    0.03267   6.707 2.15e-11 ***
body.weight   0.10341    0.01437   7.198 6.76e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 42.59 on 6734 degrees of freedom
Multiple R-squared:  0.09077,   Adjusted R-squared:  0.09037 
F-statistic: 224.1 on 3 and 6734 DF,  p-value: < 2.2e-16

# Create a copy of nhanes for imputation
#imputed_nhanes <- nhanes

# Identify rows with missing cholestrol values
missing_cholestrol <- is.na(imputed_nhanes$cholestrol)

# Use the model to predict the missing cholestrol values
imputed_nhanes$cholestrol[missing_cholestrol] <- predict(model_cholestrol, 
                                                          newdata = imputed_nhanes[missing_cholestrol, c("age", "avg.sbp", "body.weight")])

All explanatory features are significant.

# Get residuals from the model for complete_nhanes
residuals_cholestrol <- residuals(model_cholestrol)

# Predicted values for complete_nhanes
predicted_cholestrol <- predict(model_cholestrol, newdata = complete_nhanes)

# Create a data frame for residuals and predicted values
residuals_data <- data.frame(
  predicted_cholestrol = predicted_cholestrol,
  residuals_cholestrol = residuals_cholestrol
)

# Create the residual plot
library(ggplot2)
ggplot(residuals_data, aes(x = predicted_cholestrol, y = residuals_cholestrol)) +
  geom_point(color = "blue") +  # Plot the residuals
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +  # Add a horizontal line at 0
  labs(title = "Residual Plot for Cholestrol Prediction",
       x = "Predicted Cholestrol",
       y = "Residuals") +
  theme_minimal()

# Create the overlay density plot
ggplot() +
  geom_density(data = nhanes, aes(x = cholestrol), fill = "blue", alpha = 0.2, na.rm = TRUE) +  # Before imputation
  geom_density(data = imputed_nhanes, aes(x = cholestrol), fill = "red", alpha = 0.2, na.rm = TRUE) +  # After imputation
  labs(title = "Density Plot of Cholestrol: Before and After Imputation",
       x = "Cholestrol", y = "Density") +
  theme_minimal() +
  scale_fill_manual(name = "Cholestrol", values = c("blue" = "blue", "red" = "red"), 
                    labels = c("Before Imputation", "After Imputation"))

The distribution of cholesterol did not change significantly after imputation. As well, the residual plot does not display any patterns. In turn, the model fit the data well and the assumption of constant variance of residuals is satisfied.

3.3 Feature Engineer

In turn a feature which calculates BMI is created to determine observations are classified as obese. According to the CDC (https://www.cdc.gov/nccdphp/dnpao/data-trends-maps/help/npao_dtm/definitions.html), Adult obesity is defined as body mass index (BMI) ≥ 30.0. Another feature named obese, based on BMI, is created, where 0 = no and 1 = yes.

# Create BMI feature
imputed_nhanes$BMI <- (imputed_nhanes$body.weight * 703) / (imputed_nhanes$height * imputed_nhanes$height)
# Round the BMI to the nearest ones place
imputed_nhanes$BMI <- round(imputed_nhanes$BMI, 0)

# Create the obese feature based on BMI
imputed_nhanes$obese <- ifelse(imputed_nhanes$BMI >= 30, 1, 0)

par(mfrow = c(1, 2))

# Histogram for BMI

hist(imputed_nhanes$BMI, 
     main = "Histogram of BMI", 
     xlab = "BMI", 
     col = "steelblue", 
     border = "black", 
     breaks = 20)

# bar chart for obese
nhanes$obese <- factor(imputed_nhanes$obese, levels = c(0, 1), labels = c("No", "Yes"))

obese_count <- table(imputed_nhanes$obese)

barplot(height = obese_count, col = "steelblue", main = "Obesity", xlab = "Obesity", ylab = "Count")

# Set up a 2x2 plotting grid
par(mfrow = c(2, 2))


# Create a contingency table
obese_hbp_table <- table(imputed_nhanes$obese, imputed_nhanes$hbp)

# Generate mosaic plot
mosaicplot(obese_hbp_table, 
           col = c("gold", "royalblue"), 
           main = "Obesity and High Blood Pressure",
           xlab = "Obesity Status", 
           ylab = "HBP Status",
           border = "black")

# Create a contingency table
obese_sex_table <- table(imputed_nhanes$sex, imputed_nhanes$obese)

# Generate mosaic plot
mosaicplot(obese_sex_table, 
           col = c("gold", "royalblue"), 
           main = "Sex and Obesity",
           xlab = "Obesity Status", 
           ylab = "Sex",
           border = "black")


# Create a contingency table
obese_race_table <- table(imputed_nhanes$race, imputed_nhanes$obese)

# Generate mosaic plot
mosaicplot(obese_race_table, 
           col = c("gold", "royalblue"), 
           main = "Race and Obesity",
           xlab = "Obesity Status", 
           ylab = "Race",
           border = "black")


# Create a contingency table
obese_smoke_table <- table(imputed_nhanes$current.smk, imputed_nhanes$obese)

# Generate mosaic plot
mosaicplot(obese_smoke_table, 
           col = c("gold", "royalblue"), 
           main = "Current Smoking Status and Obesity",
           xlab = "Obesity Status", 
           ylab = "Current Smoking Status",
           border = "black")

3.4 Multiple Imputation

Multiple Imputation by Chained Equations (MICE) is a powerful technique used to handle missing data by creating multiple plausible imputations for each missing value based on relationships in the observed data. For body weight, cholesterol, and race, MICE can be applied to generate accurate imputations, accounting for correlations between these variables and other features in the dataset. This method allows for more reliable statistical analyses by reducing bias and preserving variability in the imputed values, making it especially useful for handling missing data in incoming new data sets.

nhanes1$race <- as.factor(nhanes1$race)  
nhanes1$cholestrol <- as.numeric(nhanes1$cholestrol)
nhanes1$body.weight <- as.numeric(nhanes1$body.weight)

nhanes2 <- nhanes1


# 1 initialization

init <- mice(nhanes2, maxit=0)
print(init$method)

          X         obs         psu     stratum stat.weight         age 
         ""          ""          ""          ""          ""          "" 
        sex        race body.weight      height     avg.sbp     avg.dbp 
         ""   "polyreg"       "pmm"          ""          ""          "" 
   past.smk current.smk  cholestrol         hbp 
         ""          ""       "pmm"          "" 

# 2 imputation models

imp <- mice(nhanes2, method = c("", "", "", "", "", "", "", "polyreg", "pmm", "", "", "", "", "", "pmm", ""), 
            maxit=10,
            m=1,
            seed = 123,
            print=F)

complete(imp, action = 1)[1:10,]

    X obs psu stratum stat.weight age sex race body.weight height avg.sbp
1   3   9   2      43    19451.83  48   0    1       149.7   61.8     131
2   5  11   2      40     1245.52  48   1    1       155.3   66.2     120
3   6  19   1      35     3860.97  44   1    2       189.6   70.2     133
4   7  34   1      13     5031.70  42   0    2       125.8   62.6     100
5  10  48   1      24    26919.29  56   0    1       239.9   67.6     128
6  12  51   1      28     2730.56  44   1    1       317.9   71.1     130
7  15  55   2      44      804.03  48   1    1       245.6   68.0     155
8  19  70   1       9     2629.39  63   1    2       202.9   67.8     137
9  20  71   1      43     1147.32  37   0    1       151.7   61.7     128
10 22  73   1      29     3991.81  42   1    3       205.0   68.4     148
   avg.dbp past.smk current.smk cholestrol hbp
1       73        1           2        236   0
2       70        1           2        260   0
3       85        1           1        187   0
4       67        1           1        216   0
5       73        1           2        156   0
6       86        1           1        162   0
7       91        1           1        212   1
8       68        1           1        186   0
9       70        1           1        209   0
10      83        1           1        267   1

# 3 multiple (iterative) imputations

imp5 <- mice(nhanes2, method = "pmm", m = 5, maxit = 10, seed = 123, print=F)
plot(imp5)

Multiple Imputation by Chained Equations (MICE) will be applied to handle missing values in the body weight variable. To assess the quality of the imputations, the Mean Squared Error (MSE) will be calculated and compared the MSE of the single imputation conducted for body weight earlier in this analysis, providing a measure of the imputation accuracy.

model5_bodyweight <- with(imp5, lm(body.weight ~ height + avg.dbp + sex))
summary.stats = summary(model5_bodyweight)


summary(pool(model5_bodyweight))

         term     estimate  std.error  statistic        df       p.value
1 (Intercept) -205.7045484 8.74935482 -23.510825 2473.3145 1.747163e-110
2      height    4.8466598 0.13536721  35.803796 1634.9416 8.934110e-208
3     avg.dbp    0.7523449 0.03941721  19.086712  157.8342  7.501512e-43
4         sex   -6.6078962 1.06330598  -6.214482  276.0479  1.892089e-09

beta = summary.stats$estimate[seq(1,15,by=3)]  # explicit vector: c(1,4,7,10,13)
beta.var = (summary.stats$std.error[seq(1,15,by=3)])^2 
Q = mean(beta)
U = mean(beta.var)
B = var(beta)
T = U + (6/5)*B
pool.se = sqrt(T)
cbind(pool.se.intercept = pool.se)

     pool.se.intercept
[1,]          123.2177

The model can be used to predict body weight for incoming new data.

The MSE for single body weight imputation used earlier in this analysis is MSE = (RSE)^2 = (33.52)^2 = 1123.6 which is substantially higher than MSE when using MICE. In turn, MICE is recommended in this analysis.

Multiple Imputation by Chained Equations (MICE) will be applied to handle missing values in the cholesterol variable. To assess the quality of the imputations, the Mean Squared Error (MSE) will be calculated and compared the MSE of the single imputation conducted for cholesterol earlier in this analysis, providing a measure of the imputation accuracy.

model5_chol <- with(imp5, lm(cholestrol ~ age + avg.sbp))
summary.stats = summary(model5_chol)


summary(pool(model5_chol))

         term    estimate  std.error statistic        df      p.value
1 (Intercept) 148.6912121 3.65965489 40.629845  74.24158 1.827860e-52
2         age   0.5233278 0.03407744 15.357017 236.74859 2.129193e-37
3     avg.sbp   0.2585339 0.03425598  7.547118  63.62530 2.094322e-10

beta = summary.stats$estimate[seq(1,15,by=3)]  # explicit vector: c(1,4,7,10,13)
beta.var = (summary.stats$std.error[seq(1,15,by=3)])^2 
Q = mean(beta)
U = mean(beta.var)
B = var(beta)
T = U + (6/5)*B
pool.se = sqrt(T)
cbind(pool.se.intercept = pool.se)

     pool.se.intercept
[1,]          3.659655

The model can be used to predict cholesterol for incoming new data.

The Mean Squared Error (MSE) from the Residual Standard Error (RSE) given in the model output from the cholesterol single imputation method used early in this analysis is MSE = (42.59)^2 = 1813.9. Thus, the Mean Squared Error (MSE) = 1813.9.

The MSE using MICE is 3.66 which is substantially less. Therefore, MICE is recommended in the analysis of this data set.

For the missing values in the race variable, Multiple Imputation by Chained Equations (MICE) will be applied using multinomial logistic regression. This method will account for the categorical nature of the race variable, allowing for appropriate imputation based on the relationships with other features in the dataset.

model_race <- with(imp5, multinom(race ~ avg.dbp + avg.sbp + cholestrol + hbp))

# weights:  18 (10 variable)
initial  value 8708.699612 
iter  10 value 5523.235718
final  value 5495.930272 
converged
# weights:  18 (10 variable)
initial  value 8708.699612 
iter  10 value 5555.595871
final  value 5522.606318 
converged
# weights:  18 (10 variable)
initial  value 8708.699612 
iter  10 value 5555.946819
final  value 5529.649391 
converged
# weights:  18 (10 variable)
initial  value 8708.699612 
iter  10 value 5560.173296
final  value 5529.439901 
converged
# weights:  18 (10 variable)
initial  value 8708.699612 
iter  10 value 5542.550592
final  value 5510.127472 
converged

# Pool the results from the imputed datasets
pooled_race <- pool(model_race)

# Get a summary of the pooled results
summary(pooled_race)

          term     estimate    std.error  statistic         df      p.value
1  (Intercept) -1.277932787 0.3325448882 -3.8428881   49.56081 3.468517e-04
2      avg.dbp  0.030728414 0.0030785469  9.9814670  343.89941 9.007490e-21
3      avg.sbp -0.009072042 0.0025775079 -3.5196952  252.21602 5.124609e-04
4   cholestrol -0.004119817 0.0007738102 -5.3240668   28.44201 1.091924e-05
5          hbp  0.148914062 0.1127433356  1.3208236  101.89054 1.895191e-01
6  (Intercept) -2.575346290 0.9149529588 -2.8147308   40.55536 7.498473e-03
7      avg.dbp  0.019877690 0.0080743228  2.4618399 3279.39165 1.387374e-02
8      avg.sbp -0.007708776 0.0071499482 -1.0781583  208.24247 2.822104e-01
9   cholestrol -0.005554712 0.0020814684 -2.6686508   36.04555 1.134283e-02
10         hbp -0.151093362 0.3197542727 -0.4725296  109.96948 6.374852e-01

The model can be used to predict race for incoming new data.

4 FEATURE TRANSFORMING

Based on the distribution patterns observed in the NHANES data set, several feature engineering tasks can be performed to enhance the data’s usability for modeling. First, for age, which is roughly uniform, no major transformation is needed, though we may consider binning it into age groups if it helps with interpretability. For weight and systolic blood pressure (BP), which are both skewed right, a log transformation could be applied to reduce skewness and make these variables more normally distributed, improving model performance. Similarly, for cholesterol, which is slightly skewed right, a log transformation may also be useful. For diastolic BP, which is roughly symmetric, no transformation is necessary, but checking for outliers could be helpful. These transformations aim to stabilize variance, reduce the impact of extreme values, and create features that are more suitable for modeling with techniques like linear regression or other statistical methods.

#Log transformation to reduce skew

imputed_nhanes$log_body.weight <- log(imputed_nhanes$body.weight)
imputed_nhanes$log_avg.sbp <- log(imputed_nhanes$avg.sbp)
imputed_nhanes$log_cholestrol <- log(imputed_nhanes$cholestrol)
imputed_nhanes$log_BMI <- log(imputed_nhanes$BMI)




# Set up the plotting area to show 2x2 grid of histograms
par(mfrow = c(2, 2))  # This will create a 2x2 grid of plots

# Histogram for Log-Cholesterol
hist(imputed_nhanes$log_cholestrol, 
     main = "Histogram of Log(Cholesterol)", 
     xlab = "Log(Cholesterol)", 
     col = "steelblue", 
     border = "black", 
     breaks = 20)  
    
# Histogram for Log-Body Weight
hist(imputed_nhanes$log_body.weight, 
     main = "Histogram of Log(Body Weight)", 
     xlab = "Log(Body Weight)", 
     col = "steelblue", 
     border = "black", 
     breaks = 20) 


# Histogram for Log-Systolic BP
hist(imputed_nhanes$log_avg.sbp, 
     main = "Histogram of Log(Systolic BP)", 
     xlab = "Log(Systolic BP)", 
     col = "steelblue", 
     border = "black", 
     breaks = 20)
   

# Histogram for Log-BMI
hist(imputed_nhanes$log_BMI, 
     main = "Histogram of Log(BMI)", 
     xlab = "Log(BMI)", 
     col = "steelblue", 
     border = "black", 
     breaks = 20)

The distributions of these previously skewed features are more symmetric after log transformations.

5 STANDARDIZING

Standardizing features ensures that variables with different scales, such as age, BMI, height, and cholesterol, contribute equally to the analysis. Without standardization, variables with larger ranges could dominate model performance, leading to biased results. By standardizing these features, we improve model accuracy, enable better comparison between variables, and enhance the efficiency of machine learning algorithms, particularly in models that rely on distance-based or gradient-based methods.

# Specify the columns to be standardized
columns_to_standardize <- c("age", "BMI", "height", "body.weight", "cholestrol", "avg.dbp", "avg.sbp")

# Standardize the selected columns
imputed_nhanes[columns_to_standardize] <- scale(imputed_nhanes[columns_to_standardize])

6 FEATURE SELECTION

Feature selection is a critical step in model building, as it helps improve model performance by identifying the most relevant predictors and reducing overfitting. In this analysis, both model-based selection and regularization-based selection techniques will be used to assess and retain the most influential features for accurate prediction, ensuring a more efficient and interpretable model.

A full model is created to predict cholesterol.

#sapply(imputed_nhanes, function(x) length(unique(x)))

#head(imputed_nhanes, 25)


#head(imputed_nhanes)


model_full <- lm(cholestrol ~ age + race + BMI + height + body.weight + 
                  hbp + avg.dbp + avg.sbp + obese, 
                  data = imputed_nhanes)

# Summary of the full model
summary(model_full)

Call:
lm(formula = cholestrol ~ age + race + BMI + height + body.weight + 
    hbp + avg.dbp + avg.sbp + obese, data = imputed_nhanes)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7355 -0.5662 -0.0433  0.4439 11.9436 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.04277    0.01748   2.447 0.014423 *  
age          0.24974    0.01325  18.847  < 2e-16 ***
raceBlack   -0.06770    0.02435  -2.781 0.005436 ** 
raceOther   -0.20562    0.06724  -3.058 0.002237 ** 
BMI         -0.12183    0.08249  -1.477 0.139735    
height      -0.18964    0.04910  -3.862 0.000113 ***
body.weight  0.26560    0.09330   2.847 0.004426 ** 
hbp         -0.06604    0.04197  -1.574 0.115618    
avg.dbp      0.11648    0.01327   8.774  < 2e-16 ***
avg.sbp      0.05361    0.02064   2.597 0.009419 ** 
obese       -0.02242    0.03806  -0.589 0.555846    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9366 on 7916 degrees of freedom
Multiple R-squared:  0.1238,    Adjusted R-squared:  0.1227 
F-statistic: 111.9 on 10 and 7916 DF,  p-value: < 2.2e-16

Stepwise model selection is conducted to retain significant predictors. Selection processes based on AIC and BIC are conducted.

# Stepwise model selection based on AIC
model_stepwise_AIC <- step(model_full, direction = "both", trace = 0)

# Stepwise model selection based on BIC
model_stepwise_BIC <- step(model_full, direction = "both", k = log(nrow(imputed_nhanes)), trace = 0)

# Display the summary of the selected model
summary(model_stepwise_AIC)

Call:
lm(formula = cholestrol ~ age + race + BMI + height + body.weight + 
    hbp + avg.dbp + avg.sbp, data = imputed_nhanes)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7302 -0.5629 -0.0437  0.4447 11.9404 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.03788    0.01538   2.462 0.013820 *  
age          0.24998    0.01324  18.875  < 2e-16 ***
raceBlack   -0.06795    0.02434  -2.792 0.005257 ** 
raceOther   -0.20591    0.06724  -3.062 0.002203 ** 
BMI         -0.12726    0.08197  -1.553 0.120563    
height      -0.18849    0.04906  -3.842 0.000123 ***
body.weight  0.26375    0.09324   2.829 0.004685 ** 
hbp         -0.06660    0.04195  -1.588 0.112433    
avg.dbp      0.11652    0.01327   8.778  < 2e-16 ***
avg.sbp      0.05375    0.02064   2.604 0.009223 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9366 on 7917 degrees of freedom
Multiple R-squared:  0.1238,    Adjusted R-squared:  0.1228 
F-statistic: 124.3 on 9 and 7917 DF,  p-value: < 2.2e-16

summary(model_stepwise_BIC)

Call:
lm(formula = cholestrol ~ age + height + body.weight + avg.dbp, 
    data = imputed_nhanes)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.6387 -0.5588 -0.0418  0.4475 11.9160 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -5.892e-16  1.053e-02   0.000        1    
age          2.714e-01  1.067e-02  25.433   <2e-16 ***
height      -1.164e-01  1.208e-02  -9.640   <2e-16 ***
body.weight  1.223e-01  1.235e-02   9.908   <2e-16 ***
avg.dbp      1.289e-01  1.119e-02  11.519   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9378 on 7922 degrees of freedom
Multiple R-squared:  0.121, Adjusted R-squared:  0.1206 
F-statistic: 272.7 on 4 and 7922 DF,  p-value: < 2.2e-16

# Check R-squared value
summary(model_stepwise_AIC)$r.squared

[1] 0.123777

summary(model_stepwise_BIC)$r.squared

[1] 0.1210105

Both AIC and BIC stepwise model selection processes have similar R squared values. The BIC stepwise model is parsimonious. This model selects age, height, body weight, and average diastolic blood pressure as predictors for cholesterol.

Regularization-based selection techniques are useful for improving model performance and preventing overfitting, especially when you have many predictors. Two common regularization methods are Ridge regression and Lasso regression.

# Extract the response and predictors
y <- imputed_nhanes$cholestrol  # Response variable
X <- as.matrix(imputed_nhanes[, c("age", "BMI", "height", "body.weight", "avg.dbp", "avg.sbp")])

# Optional: Scale predictors (recommended for Ridge and Lasso)
X_scaled <- scale(X)

ridge_model <- glmnet(X_scaled, y, alpha = 0)  # alpha = 0 for Ridge

lasso_model <- glmnet(X_scaled, y, alpha = 1)  # alpha = 1 for Lasso

# Ridge cross-validation
cv_ridge <- cv.glmnet(X_scaled, y, alpha = 0)
plot(cv_ridge)

best_lambda_ridge <- cv_ridge$lambda.min

# Lasso cross-validation
cv_lasso <- cv.glmnet(X_scaled, y, alpha = 1)
plot(cv_lasso)

best_lambda_lasso <- cv_lasso$lambda.min


# Ridge coefficients
ridge_coefficients <- coef(cv_ridge, s = "lambda.min")
print(ridge_coefficients)

7 x 1 sparse Matrix of class "dgCMatrix"
                       s1
(Intercept)  4.291996e-15
age          2.453805e-01
BMI          4.128314e-02
height      -8.611899e-02
body.weight  7.098814e-02
avg.dbp      1.085174e-01
avg.sbp      3.889874e-02

# Lasso coefficients
lasso_coefficients <- coef(cv_lasso, s = "lambda.min")
print(lasso_coefficients)

7 x 1 sparse Matrix of class "dgCMatrix"
                       s1
(Intercept)  4.290261e-15
age          2.554719e-01
BMI          .           
height      -1.128253e-01
body.weight  1.196546e-01
avg.dbp      1.139278e-01
avg.sbp      3.029412e-02

# Predictions for Ridge
pred_ridge <- predict(cv_ridge, X_scaled, s = "lambda.min")
# Predictions for Lasso
pred_lasso <- predict(cv_lasso, X_scaled, s = "lambda.min")

# Calculate RMSE or other metrics to assess model performance

# Predictions for Ridge
pred_ridge <- predict(cv_ridge, X_scaled, s = "lambda.min")

# Predictions for Lasso
pred_lasso <- predict(cv_lasso, X_scaled, s = "lambda.min")

# Calculate RMSE for Ridge
rmse_ridge <- sqrt(mean((pred_ridge - y)^2))

# Calculate RMSE for Lasso
rmse_lasso <- sqrt(mean((pred_lasso - y)^2))

# Print RMSE values
cat("Ridge RMSE:", rmse_ridge, "\n")

Ridge RMSE: 0.9374275 

cat("Lasso RMSE:", rmse_lasso, "\n")

Lasso RMSE: 0.9372577 

# Calculate R² for Ridge
rss_ridge <- sum((pred_ridge - y)^2)  # Residual Sum of Squares
tss_ridge <- sum((y - mean(y))^2)    # Total Sum of Squares
r2_ridge <- 1 - rss_ridge / tss_ridge

# Calculate R² for Lasso
rss_lasso <- sum((pred_lasso - y)^2)
tss_lasso <- sum((y - mean(y))^2)
r2_lasso <- 1 - rss_lasso / tss_lasso

# Print R² values
cat("Ridge R²:", r2_ridge, "\n")

Ridge R²: 0.1211188 

cat("Lasso R²:", r2_lasso, "\n")

Lasso R²: 0.1214371 

The Ridge and Lasso method both had similar R squared values. Therefore, both models are recommended.